home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 001 / pn.arc / PN.PLB next >
Text File  |  1985-08-12  |  7KB  |  168 lines
  1. {      PN.PLB #2.00 85-08-06 NORMAL PROBABILITY DISTRIBUTION FUNCTION
  2.  
  3.        V02 L00 85-08-06 conversion to Turbo Pascal by Dennis E. Hamilton,
  4.                for separate testing for use as part of the Chi-Squared
  5.                Distribution procedure, Q2(x, df).  This version is actually
  6.                a clone of ERF.PLB #2.00, since the functions are so similar.
  7.  
  8.        V01 L01 79-03-22 implementated in Benton Harbor BASIC 05.01.01 by DEH;
  9.                used in deriving chi-squared distributions in testing RND().
  10.            L00 79-03-18 Texas Instruments SR52 calculator program derived
  11.                by Dennis E. Hamilton to verify basic recurrence methods and
  12.                checking against values in standard tables.       }
  13.  
  14.  
  15. function
  16.  
  17.    PN(x: real {value of normally-distributed random variable} )
  18.  
  19.       :real {approximate Normal Probability Distribution Function at x};
  20.  
  21.       {Probability that, for Y a normally-distributed random variable
  22.        with mean 0 and standard deviation 1.0, Y will take on  a value
  23.        not greater than x.  100*PN(x) is the percentile of the entire
  24.        Y population having value not greater than x.  Similarly, QN(x),
  25.        defined to be 1-PN(x), is the probability with which Y > x.       }
  26.  
  27.  
  28.    const
  29.  
  30.       sqrt2 = 1.414213562373;
  31.               {an accurate-enough value of sqrt(2) from [HMF: Table 1.1]}
  32.  
  33.       e1max = 1E-10;
  34.               {largest positive value of x for which exp(x) is exactly 1.0
  35.                to the floating-point precision of this computer.  This value
  36.                is used to avoid computation of insignificant terms.}
  37.  
  38.       e0min = 81;
  39.               {smallest nonzero value of x at which exp(-x) underflows to
  40.                0.0 within the computer's floating-point precision.  This value
  41.                is used to steer around avoidable arithmetic overflows.}
  42.  
  43.  
  44.    var y: real {intermediate result};
  45.  
  46.    function H1(x: real)
  47.  
  48.                :real {approximation of 1-erf(abs(x)) where erf(x) is
  49.                       the standard Gaussian error integral [HMF: 7.1.1]};
  50.  
  51.       const   a6 = 0.0000430638;      a5 = 0.0002765672;
  52.               a4 = 0.0001520143;      a3 = 0.0092705272;
  53.               a2 = 0.0422820123;      a1 = 0.0705230784;
  54.  
  55.               xmax = 9 {value of sqrt(e0min)};
  56.  
  57.       begin {Computation to over 6 decimal-place precision based on Hastings
  58.              approximation 7.1.28 in [Abramowitz, Milton., Stegun, Irene A.
  59.              (eds.)  HANDBOOK OF MATHEMATICAL FUNCTIONS.  National Bureau of
  60.              Standards Applied Mathematics Series Publication #55.  Dover
  61.              Publications (New York: 1965)].  For more about this function,
  62.              consult the description of companion function erf(x) in file
  63.              ERF.PLB #2.00.    }
  64.  
  65.       x := abs(x);
  66.       if x > xmax
  67.       then H1 := 0
  68.       else begin
  69.            x := (((((a6*x + a5)*x + a4)*x + a3)*x + a2)*x + a1)*x + 1.0;
  70.            H1 := sqr(sqr(sqr(sqr(1.0/x))));
  71.            end;
  72.       end {H1};
  73.  
  74.  
  75.    BEGIN {PN};
  76.  
  77.          {Approximation of the Normal (Gaussian) Probability Distribution
  78.           Function based on standard relationships [HMF: section 26.2] and
  79.           availability of approximation H1(x) = 1 - erf(abs(x)).
  80.  
  81.           To also employ erf(x) and H1(x) independent of PN(x) computation,
  82.           remove H1(x) from here, making it global and available for separate
  83.           usage by PN, erf, and any other functions that depend on it.   }
  84.  
  85.  
  86.    y := H1(x/sqrt2)/2.0;    {PN(-x) = 1 - PN(x),             [HMF 26.2.5-2.6]}
  87.    if x < 0                 { PN(x) = (1+erf(x/sqrt2))/2,    [HMF 7.1.22]    }
  88.    then PN := y             {      = 1 - H1(x/sqrt2)/2,      x not negative  }
  89.    else PN := 1.0 - y;      {      = H1(x/sqrt2)/2,          x < 0           }
  90.  
  91.    END {PN};
  92.  
  93. {  The following results were obtained with a test program (PNT1.PAS) using
  94.    test data file PNT1.DAT #1.00.  PNT1 was compiled and executed using 
  95.    Borland International Turbo Pascal 3.00A for CP/M-80.  This version 
  96.    maintains a 40-bit floating-point significand and decimal output 
  97.    conversion rounded to 11 significant digits.  The error terms represent
  98.    the absolute difference between computed values and those given in
  99.    known tables.  In checking behavior at important extreme points, the
  100.    results of tests with ERF.PLB #2.00 were used to obtain trial values.
  101.  
  102.   
  103. PNT1> #2.00 85-08-06 NORMAL PROBABILITY DISTRIBUTION FUNCTION
  104.  
  105.             x                  PN(x)           abs(error)
  106.  
  107.        -1.273000000E+01   0.0000000000E+00      0.00E+00
  108.        -1.272000000E+01   1.8175366526E-28      1.82E-28
  109.  
  110.        -9.013270000       4.5747132329E-18      4.47E-18
  111.        -7.348800000       2.9209847295E-13      1.92E-13
  112.  
  113.        -6.706020000       1.7491994026E-11      7.49E-12
  114.        -5.997810000       1.2641232113E-09      2.64E-10
  115.        -5.199340000       0.00000010694         6.94E-09
  116.        -4.753420000       0.00000102818         2.82E-08
  117.        -4.264890000       0.00001008388         8.39E-08
  118.        -3.719020000       0.00010012594         1.26E-07
  119.  
  120.        The values computed for PN(x),  with x very small,  are
  121.        compared with tabulated values of PN(-x) = QN(x) in The
  122.        Handbook of  Mathematical Functions  [HMF: Table 26.6].
  123.        Errors within 1E-5 are usually quite acceptable.   Note
  124.        the very high relative errors out beyond -6.0, however.
  125.  
  126.             z                  PN(x)           abs(error)
  127.  
  128.        -0.140000000       0.44433009345         9.82E-08
  129.        -0.020000000       0.49202174538         5.91E-08
  130.        -0.002510000       0.49899866454         1.34E-06
  131.  
  132.        -3.647656264E-11   0.49999999999         1.32E-11
  133.        -3.647514842E-11   0.50000000000         0.00E+00
  134.         0.000000000E+00   0.50000000000         0.00E+00
  135.         3.647514842E-11   0.50000000000         0.00E+00
  136.         3.647656264E-11   0.50000000001         1.36E-11
  137.  
  138.         0.002510000       0.50100133546         1.34E-06
  139.         0.020000000       0.50797825462         5.91E-08
  140.         0.140000000       0.55566990655         9.83E-08
  141.  
  142.        Midrange values, near .5, check on the smoothness of
  143.        the crossover from negative to positive values of x.
  144.        Test values are from [HMF:Table 26.1] and ERFT1.PAS.
  145.  
  146.             x                  PN(x)           abs(error)
  147.  
  148.         2.326350000       0.99000003622         3.62E-08
  149.         3.090230000       0.99900004431         4.43E-08
  150.         3.719020000       0.99989987406         1.26E-07
  151.         4.264890000       0.99998991612         8.39E-08
  152.         4.753420000       0.99999897182         2.82E-08
  153.         5.199340000       0.99999989307         6.93E-09
  154.         5.612000000       0.99999998857         1.43E-09
  155.         5.997810000       0.99999999874         2.62E-10
  156.         6.361340000       0.99999999986         4.37E-11
  157.         6.706020000       0.99999999998         6.37E-12
  158.         7.034480000       1.00000000000         0.00E+00
  159.  
  160.        Near the tail, as PN(x) approaches 1.0,  values of 1-QN(x)
  161.        from [HMF: Table 26.6] are used to see how well PN(x) does
  162.        despite adjustments by 1.0 within the PN procedure.
  163.  
  164. PNT1> End of test.
  165.                                 }
  166.             (* end of PN.PLB *)
  167.  
  168.